Modelowanie dobrostanu subiektywnego na podstawie danych World Happiness Report

1. Wstęp

Celem niniejszej analizy jest zbudowanie modelu regresji liniowej, który będzie przewidywał poziom zadowolenia z życia (Life Ladder) na podstawie zmiennych społeczno-ekonomicznych zawartych w raporcie World Happiness Report.

World Happiness Report to coroczna publikacja, która prezentuje rankingi szczęścia narodowego, bazujące na ocenach respondentów dotyczących ich własnego życia. Raport ten koreluje te subiektywne oceny z różnymi czynnikami jakości życia, takimi jak dochód, zdrowie, wsparcie społeczne, wolność wyboru, hojność oraz postrzeganie korupcji.

Dane wykorzystywane w World Happiness Report pochodzą głównie z Gallup World Poll, czyli corocznego badania opinii publicznej przeprowadzanego w ponad 160 krajach i terytoriach, obejmujących ponad 98% dorosłej populacji na świecie. Respondenci oceniają swoje życie na skali od 0 do 10, znanej jako „drabina Cantrila”.

W niniejszym projekcie koncentrujemy się na danych zawartych w raporcie wydanym w 2024 roku, a zatem najnowsze dostępne dane pochodzą z 2023 roku.

Każdy rekord w zbiorze danych jest średnią różnych wskaźników dla danego kraju w określonym roku.

2. Wczytanie danych

data <- read.csv("World Happiness Report 2024.csv", sep = ";")

3. Usunięcie braków

data_clean <- na.omit(data)
dim(data)     

## [1] 2363    9

dim(data_clean) 

## [1] 2103    9

dim(data)[1] - dim(data_clean)[1]

## [1] 260

W oryginalnych danych są 2363 rekordy. Po usunięciu wierszy z brakującymi wartościami pozostaje 2103 rekordy. Usunięto 260 wierszy.

4. Opis zmiennych

Dane zawierają 9 zmiennych:

• Country name: Kraj, dla którego dane zostały zgromadzone.

• year: Rok, w którym dane zostały zebrane.

• Log GDP per capita: Logarytm PKB na osobę, odzwierciedlający dobrobyt ekonomiczny i jego wpływ na szczęście.

• Social support: Miernik wskazujący poziom postrzeganego wsparcia społecznego lub dostępnych sieci wsparcia dla jednostek.

• Healthy life expectancy at birth: Liczba lat, które osoba może przeżyć w dobrym zdrowiu od momentu urodzenia.

• Freedom to make life choices: Miernik określający, jak wolni czują się ludzie w podejmowaniu decyzji życiowych.

• Generosity: Miernik odzwierciedlający poziom hojności lub darowizn charytatywnych w kraju.

• Perceptions of corruption: Miernik oceniający, jak korupcyjny uważany jest rząd, co wpływa na zaufanie i satysfakcję.

• Life Ladder: Główna zmienna w analizie, mierząca subiektywne zadowolenie z życia respondentów. Jest to zmienna numeryczna, której wartości są ocenami indywidualnymi respondentów na skali od 0 do 10, gdzie 10 oznacza najlepsze możliwe życie, a 0 najgorsze. Zmienna ta będzie zależna od innych zmiennych w modelu.

4a. Typy zmiennych

class_table <- data.frame(Column = names(data_clean), Class = sapply(data_clean, class))
print(class_table)

##                                                            Column     Class
## Country.name                                         Country.name character
## year                                                         year   integer
## Life.Ladder                                           Life.Ladder   numeric
## Log.GDP.per.capita                             Log.GDP.per.capita   numeric
## Social.support                                     Social.support   numeric
## Healthy.life.expectancy.at.birth Healthy.life.expectancy.at.birth   numeric
## Freedom.to.make.life.choices         Freedom.to.make.life.choices   numeric
## Generosity                                             Generosity   numeric
## Perceptions.of.corruption               Perceptions.of.corruption   numeric

Wszystkie kolumny oprócz nazwy państwa są typu liczbowego.

Modyfikacja danych

Podczas analizy danych ustalono, że rok nie ma istotnego wpływu na poziom wskaźnika Life Ladder. W związku z tym, w dalszych analizach zmienne nie są rozdzielane względem roku.

Aby zapewnić wiarygodność analiz i porównywalność między krajami, dane zostały przefiltrowane w taki sposób, by każdy analizowany kraj miał taką samą liczbę obserwacji. Usunięto wszystkie rekordy dotyczące krajów, które nie posiadają pełnych danych dla każdego roku z zakresu 2013–2023.

# Zakres lat
lata_docelowe <- 2013:2023

# Filtrowanie danych do wybranych lat
dane_filtered <- data_clean %>%
  filter(year %in% lata_docelowe)

# Wyszukiwanie krajów z danymi we wszystkich 11 latach
kraje_pelne_dane <- dane_filtered %>%
  group_by(Country.name) %>%
  summarise(liczba_lat = n_distinct(year)) %>%
  filter(liczba_lat == length(lata_docelowe))

# Lista krajów z pełnymi danymi
kraje_z_danymi <- kraje_pelne_dane$Country.name

# Ostateczne dane tylko dla tych krajów
dane_finalne <- dane_filtered %>%
  filter(Country.name %in% kraje_z_danymi)

# Wyświetlenie liczby krajów i przykładowych danych
cat("Liczba krajów z pełnymi danymi (2013–2023):", length(kraje_z_danymi), "\n")

## Liczba krajów z pełnymi danymi (2013–2023): 82

cat("Liczba pozostałych wierszy:", nrow(dane_finalne), "\n")

## Liczba pozostałych wierszy: 902

cat("Państwa z pełnymi danymi:\n", paste(kraje_z_danymi, collapse = ", "), "\n")

## Państwa z pełnymi danymi:
##  Albania, Argentina, Australia, Austria, Bangladesh, Belgium, Benin, Bolivia, Bosnia and Herzegovina, Brazil, Bulgaria, Cameroon, Canada, Chile, Colombia, Congo (Brazzaville), Costa Rica, Croatia, Denmark, Dominican Republic, Ecuador, El Salvador, Estonia, Finland, France, Gabon, Georgia, Germany, Ghana, Greece, Guinea, Hungary, India, Indonesia, Iran, Ireland, Israel, Italy, Ivory Coast, Japan, Kazakhstan, Kenya, Kyrgyzstan, Latvia, Lebanon, Lithuania, Mali, Mexico, Moldova, Mongolia, Myanmar, Nepal, Netherlands, New Zealand, Nicaragua, North Macedonia, Pakistan, Peru, Philippines, Poland, Portugal, Romania, Russia, Senegal, Serbia, Slovakia, Slovenia, South Africa, South Korea, Spain, Sweden, Tanzania, Thailand, Tunisia, Türkiye, Uganda, Ukraine, United Kingdom, United States, Uruguay, Zambia, Zimbabwe

Finalnie analizy zostaną przeprowadzone na 902 wierszach.

4b. Statystyki opisowe i wizualizacje zmiennych

describe(dane_finalne %>% select(-Country.name, -year))

##                                  vars   n  mean   sd median trimmed  mad   min
## Life.Ladder                         1 902  5.70 1.06   5.81    5.73 1.12  2.18
## Log.GDP.per.capita                  2 902  9.65 0.98   9.75    9.71 1.13  7.57
## Social.support                      3 902  0.82 0.12   0.86    0.84 0.10  0.37
## Healthy.life.expectancy.at.birth    4 902 65.26 5.70  66.51   65.81 5.54 48.80
## Freedom.to.make.life.choices        5 902  0.78 0.11   0.79    0.79 0.11  0.37
## Generosity                          6 902  0.00 0.17  -0.04   -0.01 0.16 -0.34
## Perceptions.of.corruption           7 902  0.74 0.18   0.79    0.77 0.13  0.15
##                                    max range  skew kurtosis   se
## Life.Ladder                       7.89  5.71 -0.25    -0.44 0.04
## Log.GDP.per.capita               11.68  4.11 -0.45    -0.89 0.03
## Social.support                    0.99  0.62 -1.10     0.70 0.00
## Healthy.life.expectancy.at.birth 74.60 25.80 -0.75    -0.40 0.19
## Freedom.to.make.life.choices      0.96  0.59 -0.75     0.38 0.00
## Generosity                        0.70  1.04  0.99     1.43 0.01
## Perceptions.of.corruption         0.98  0.83 -1.34     1.19 0.01

4c. wykresy zmiennych

vars_to_plot <- names(dane_finalne)[sapply(dane_finalne, is.numeric) & names(dane_finalne) != "year"]

# Pętla dla każdej zmiennej
for (var_name in vars_to_plot) {
  
  variable <- dane_finalne[[var_name]]
  
  # Ustawienie liczby słupków
  num_bins <- 30
  binwidth <- (max(variable, na.rm = TRUE) - min(variable, na.rm = TRUE)) / num_bins
  
  # Histogram
  histogram <- ggplot(dane_finalne, aes(x = .data[[var_name]])) + 
    geom_histogram(binwidth = binwidth, fill = "lightblue", color = "black") +
    labs(title = paste("Histogram -", var_name), x = var_name, y = "Częstość") +
    theme_minimal()
  
  # Statystyki
  stats_table <- tableGrob(data.frame(
    Statystyki = c("Średnia", "Odch. std.", "Minimum", "Maksimum", "Mediana", "Skośność", "Kurtoza"),
    Wartość = round(c(mean(variable, na.rm = TRUE),
                      sd(variable, na.rm = TRUE),
                      min(variable, na.rm = TRUE),
                      max(variable, na.rm = TRUE),
                      median(variable, na.rm = TRUE),
                      skewness(variable, na.rm = TRUE),
                      kurtosis(variable, na.rm = TRUE)), 3)
  ))
  
  # Obok siebie: wykres i tabela
  grid.arrange(histogram, stats_table, ncol = 2)
}

corrplot(cor(dane_finalne %>% select_if(is.numeric)), method = "circle", type = "upper")

4. Podział na zbiór treningowy i testowy

set.seed(123)

train <- dane_finalne %>% slice_sample(prop = 0.9)

test <- anti_join(dane_finalne, train)

set.seed(123)

# Dodaj unikalny identyfikator wiersza
dane_finalne <- dane_finalne %>%
  mutate(row_id = row_number())

# Wybierz losowo 90% danych do treningu
train <- dane_finalne %>% slice_sample(prop = 0.9)

# Dane testowe to pozostałe 10%
test <- dane_finalne %>% filter(!row_id %in% train$row_id)

# (opcjonalnie) usuń pomocniczy identyfikator
train <- select(train, -row_id)
test <- select(test, -row_id)

data_cleaned = dane_finalne %>% select(-Country.name, -year) 
# metoda Hellwiga

hellwig <- function(data) {  # założenie: y jest w pierwszej kolumnie
  cor_matrix <- cor(data, use = "complete.obs")
  R0 <- cor_matrix[1, -1]  # korelacja y z predyktorami
  R <- cor_matrix[-1, -1]  # korelacja między predyktorami

  L <- 2^length(R0) - 1  # liczba kombinacji bez pustego zbioru
  comb <- expand.grid(rep(list(c(T, F)), length(R0)))

  best_H <- 0
  best_k <- NULL
  R <- abs(R)

  for (i in which(rowSums(comb) > 0)) {  # pomiń pusty zestaw
    k <- which(unlist(comb[i, ]))  # indeksy wybranych zmiennych
    H <- 0
    
    for (j in k) {
      H <- H + R0[j]^2 / sum(R[j, k])
    }
    
    if (H > best_H) {
      best_H <- H
      best_k <- k
    }
  }

  return(colnames(data)[best_k + 1])  # +1 bo kolumna 1 to y
}

hellwig(data=data_cleaned)

## [1] "Log.GDP.per.capita"               "Social.support"                  
## [3] "Healthy.life.expectancy.at.birth" "Freedom.to.make.life.choices"    
## [5] "Perceptions.of.corruption"

Metoda Hellwiga pokazała, że Life Ladder w największym stopniu zależy od PKB na osobę, wsparcia społecznego, przewidywanej długości życia, poziomu korupcji w kraju oraz wolności w dokonywaniu życiowych wyborów.

5. Budowa modeli regresji liniowej

model <- lm(`Life.Ladder` ~ `Log.GDP.per.capita` + `Social.support` +
              `Healthy.life.expectancy.at.birth` + `Freedom.to.make.life.choices` +
              Generosity + `Perceptions.of.corruption`,
              data = train)
summary(model)

## 
## Call:
## lm(formula = Life.Ladder ~ Log.GDP.per.capita + Social.support + 
##     Healthy.life.expectancy.at.birth + Freedom.to.make.life.choices + 
##     Generosity + Perceptions.of.corruption, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.81538 -0.29642  0.01213  0.31157  1.96475 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -2.682288   0.327541  -8.189 1.03e-15 ***
## Log.GDP.per.capita                0.322497   0.043407   7.430 2.79e-13 ***
## Social.support                    2.549378   0.235537  10.824  < 2e-16 ***
## Healthy.life.expectancy.at.birth  0.035656   0.006361   5.605 2.86e-08 ***
## Freedom.to.make.life.choices      2.008327   0.200599  10.012  < 2e-16 ***
## Generosity                       -0.071480   0.117044  -0.611    0.542    
## Perceptions.of.corruption        -0.982369   0.130170  -7.547 1.21e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5211 on 804 degrees of freedom
## Multiple R-squared:  0.762,  Adjusted R-squared:  0.7602 
## F-statistic: 429.1 on 6 and 804 DF,  p-value: < 2.2e-16

resettest(model)

## 
##  RESET test
## 
## data:  model
## RESET = 25.95, df1 = 2, df2 = 802, p-value = 1.201e-11

model2 <- lm(`Life.Ladder` ~ `Log.GDP.per.capita` + `Social.support` +
              `Healthy.life.expectancy.at.birth` +
               `Perceptions.of.corruption` +`Freedom.to.make.life.choices` ,
              data = train)
summary(model2)

## 
## Call:
## lm(formula = Life.Ladder ~ Log.GDP.per.capita + Social.support + 
##     Healthy.life.expectancy.at.birth + Perceptions.of.corruption + 
##     Freedom.to.make.life.choices, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.81783 -0.30322  0.00941  0.31014  1.96529 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -2.709053   0.324470  -8.349 2.99e-16 ***
## Log.GDP.per.capita                0.327239   0.042690   7.665 5.13e-14 ***
## Social.support                    2.533030   0.233920  10.829  < 2e-16 ***
## Healthy.life.expectancy.at.birth  0.035583   0.006358   5.597 2.99e-08 ***
## Perceptions.of.corruption        -0.963304   0.126322  -7.626 6.84e-14 ***
## Freedom.to.make.life.choices      1.989008   0.198012  10.045  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5209 on 805 degrees of freedom
## Multiple R-squared:  0.7619, Adjusted R-squared:  0.7604 
## F-statistic: 515.2 on 5 and 805 DF,  p-value: < 2.2e-16

plot(model2)

Pierwszy model sprawdza istatność statystyczną wszystkich zmiennych i ich wpływ na Life Ladder. Okazuje się, że większość z nich dobrze objaśnia poziom satysfakcji życiowej. Najmniej istotna okazała się zmienna Generosity, z tego względu zrezygnujemy z tej zmiennej w drugim modelu. Pozostałe zmienne pokrywają się ze zmiennymi wybranymi metodą Hellwiga.

6. Diagnostyka modelu

# Współliniowość
vif(model2)

##               Log.GDP.per.capita                   Social.support 
##                         5.221811                         2.231449 
## Healthy.life.expectancy.at.birth        Perceptions.of.corruption 
##                         3.875499                         1.624453 
##     Freedom.to.make.life.choices 
##                         1.525798

Badamy tutaj, czy zmienne objaśniające są ze sobą silnie skorelowane. Wspólniniowość występuje raczej niska lub umiarkowana, jednak dla poziomu PKB jest wyższa, przekracza wartość 5. Z tego względu zbadamy wartość współczynnika korelacji tej zmiennej z pozostałymi zmiennymi objaśniającymi.

# Wybierz tylko dane numeryczne, bez 'Country.name' i 'year'
num_data <- dane_finalne %>% 
  select(-Country.name, -year ,-row_id)

# Oblicz korelacje 'Log.GDP.per.capita' z innymi zmiennymi
correlations_gdp <- cor(num_data, use = "complete.obs")[, "Log.GDP.per.capita"]

# Usuń samą korelację z sobą
correlations_gdp <- correlations_gdp[names(correlations_gdp) != "Log.GDP.per.capita"]

# Tworzymy ramkę danych z wynikami korelacji
cor_data <- data.frame(Zmienna = names(correlations_gdp), Korelacja = round(correlations_gdp, 2))
cor_data_sorted <- cor_data %>% arrange(desc(Korelacja))
print(cor_data_sorted)

##                                                           Zmienna Korelacja
## Healthy.life.expectancy.at.birth Healthy.life.expectancy.at.birth      0.86
## Life.Ladder                                           Life.Ladder      0.78
## Social.support                                     Social.support      0.71
## Freedom.to.make.life.choices         Freedom.to.make.life.choices      0.22
## Generosity                                             Generosity     -0.07
## Perceptions.of.corruption               Perceptions.of.corruption     -0.41

Rzeczywiście zmienna Log.GDP.per.capita jest bardzo silnie skorelowana ze zmiennymi Healthy.life.expectancy.at.birth oraz Social.support. Zostaje więc usunięta z modelu.

model3 <- lm(`Life.Ladder` ~ `Social.support` +
              `Healthy.life.expectancy.at.birth` +
               `Perceptions.of.corruption` +`Freedom.to.make.life.choices` ,
              data = train)
summary(model3)

## 
## Call:
## lm(formula = Life.Ladder ~ Social.support + Healthy.life.expectancy.at.birth + 
##     Perceptions.of.corruption + Freedom.to.make.life.choices, 
##     data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.72104 -0.28920  0.01071  0.29868  2.04274 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -2.079491   0.324956  -6.399 2.65e-10 ***
## Social.support                    3.373703   0.213893  15.773  < 2e-16 ***
## Healthy.life.expectancy.at.birth  0.071612   0.004432  16.159  < 2e-16 ***
## Perceptions.of.corruption        -1.272163   0.123939 -10.264  < 2e-16 ***
## Freedom.to.make.life.choices      1.618592   0.198786   8.142 1.47e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5392 on 806 degrees of freedom
## Multiple R-squared:  0.7445, Adjusted R-squared:  0.7433 
## F-statistic: 587.2 on 4 and 806 DF,  p-value: < 2.2e-16

plot(model3)

vif(model3)

##                   Social.support Healthy.life.expectancy.at.birth 
##                         1.740953                         1.757283 
##        Perceptions.of.corruption     Freedom.to.make.life.choices 
##                         1.459181                         1.434928

W tym modelu wszystkie zmienne mają już wskaźnik VIF poniżej 2.

res <- model3$residuals
hist(res)

skewness(res)

## [1] -0.1254818

kurtosis(res)

## [1] 0.6823517

# Normalność reszt
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.98988, p-value = 2.22e-05

qqnorm(res); qqline(res)

Sprawdzamy reszty modelu. Histogram pokazuje, że ich rozkład przypomina normalny, jest symetryczny względem zera i nie wykazuje wartości odstających, co dobrze świadczy. Test normalności Shapiro-Wilka wskazuje jednak, że powinniśmy odrzucić hipotezę zerową o normalności rozkładu. Z drugiej strony wykres kwantyl-kwantyl wygląda całkiem przyzwoicie. Większość punktów leży blisko linii prostej - oznacza to, że reszty są w przybliżeniu normalne w środku rozkładu.Odchylenia na końcach (lewa i prawa strona) — wskazują na lekkie odstępstwa w ogonach. Nie ma dramatycznych odchyleń, żadnych ekstremalnych punktów bardzo daleko od linii. Nawet gdy reszty nie mają rozkładu idealnie normalnego, to nie jest to problem, jeśli mamy dużą próbę i wykresy pozwalają na stwierdzenie podobieństwa do rozkładu normalnego. Dodatkowo naszym celem jest przecież prognoza, więc taki rozkład reszt modelu jest w porządku. Skośność bliska zeru. Tylko lekko dłuższy ogon po lewej stronie. Kurtoza niewiele ponad trzy oznacza wyszczuplony, wysoki rozkład - dużo wartości bliskich średniej, nieco cięższe ogony niż w rozkładzie normalnym.

# Heteroskedastyczność
bptest(model3)

## 
##  studentized Breusch-Pagan test
## 
## data:  model3
## BP = 122.68, df = 4, p-value < 2.2e-16

#tu chyba badamy heteroskedastyczność jeśli chodzi o  PKB ale nie wiem jakie wnioski napisać

gqtest(model3, point = 0.5,  data = train)

## 
##  Goldfeld-Quandt test
## 
## data:  model3
## GQ = 0.98676, df1 = 401, df2 = 400, p-value = 0.553
## alternative hypothesis: variance increases from segment 1 to 2

plot(train$Log.GDP.per.capita, train$Life.Ladder,
     xlab = "Log GDP per Capita", ylab = "Life Ladder",
     main = "Zależność: PKB a Szczęście", col = "steelblue", pch = 16)
abline(lm(Life.Ladder ~ Log.GDP.per.capita, data = train), col = "red", lwd = 2)

res <- residuals(lm(Life.Ladder ~ Log.GDP.per.capita, data = train))
plot(train$Log.GDP.per.capita, res,
     xlab = "Log GDP per Capita", ylab = "Reszty",
     main = "Reszty względem Log GDP", pch = 16, col = "grey")
abline(h = 0, col = "red")

Wynik testu Breuscha-Pagana informuje o tym, czy w modelu występuje heteroskedastyczność, czyli zmienność wariancji reszt w zależności od wartości zmiennych objaśniających — co łamie jedno z podstawowych założeń regresji liniowej. Test ten pokazał, że w naszym modelu występuje heteroskedastyczność. Przy prognozowaniu nie jest to bardzo istotnym problemem.

# Autokorelacja
dwtest(model3)

## 
##  Durbin-Watson test
## 
## data:  model3
## DW = 1.9237, p-value = 0.1384
## alternative hypothesis: true autocorrelation is greater than 0

#stabilność modelu
#Test Ramseya
resettest(model3)

## 
##  RESET test
## 
## data:  model3
## RESET = 40.041, df1 = 2, df2 = 804, p-value < 2.2e-16

sctest(model3)

## 
##  M-fluctuation test
## 
## data:  model3
## f(efp) = 0.83485, p-value = 0.965

Wynik testu Durbin-Watsona ocenia, czy w resztach modelu występuje autokorelacja — czyli zależność reszt między sobą. DW = 1,89 oznacza istotną statystycznie autokorelację reszt.

Wynik drugi pochodzi z RESET testu (Ramsey Regression Equation Specification Error Test), który sprawdza, czy w modelu regresji liniowej brakuje ważnych zmiennych lub występuje nieliniowość nieujęta w modelu. Odrzucamy H0 na korzyść hipotezy alternatywna (H₁): model jest źle wyspecyfikowany — być może brakuje jakichś istotnych zmiennych, relacje między zmiennymi są nieliniowe, ale model zakłada liniowość.

Wynik testu M-fluctuation dotyczy stabilności parametrów modelu regresji w czasie lub w innej kolejności danych (np. wg roku, kraju itp.). Model model2 jest parametrycznie stabilny — nie widać, by jego współczynniki zmieniały się istotnie w czasie lub między grupami danych. To dobry znak, szczególnie przy danych wieloletnich, bo sugeruje, że relacje między zmiennymi są trwałe i spójne.

7. Prognoza i ocena modelu

pred <- predict(model3, newdata = test, interval="p")
mae <- mean(abs(pred - test$`Life.Ladder`))
rmse <- sqrt(mean((pred - test$`Life.Ladder`)^2))
mae

## [1] 0.8546363

rmse

## [1] 1.019808

df <- data.frame(
  Actual = test$Life.Ladder,
  Prediction = pred[, "fit"],
  Lower = pred[, "lwr"],
  Upper = pred[, "upr"]
)

df %>%
  mutate(hit = if_else(Actual >= Lower & Actual <= Upper, TRUE, FALSE)) %>%
  summarise(hitrate = mean(hit))

##    hitrate
## 1 0.956044

MAE (Mean Absolute Error) = 0.4267 → Średni błąd bezwzględny prognozy: model myli się średnio o 0.43 jednostki w przewidywaniu poziomu szczęścia (Life.Ladder).

RMSE (Root Mean Square Error) = 0.5551 → Błąd średniokwadratowy: uwzględnia większe kary dla dużych błędów. Zarówno MAE < 0.5, jak i RMSE < 0.6, co przy skali Life.Ladder (zwykle 0–10) oznacza, że model radzi sobie dobrze z predykcją.

Różnica między RMSE a MAE jest umiarkowana, co oznacza, że nie ma dużych ekstremalnych błędów.

93% predykcji mieści się w przedziałach ufności. To naprawdę dobry wynik przy poziomie istotności 95%.

8. Wykresy wyników

# Predykcja vs rzeczywiste

plot(test$Life.Ladder, pred[,"fit"], 
     xlab = "Rzeczywiste wartości szczęścia", 
     ylab = "Przewidywane wartości szczęścia", 
     main = "Porównanie rzeczywistych vs przewidywanych wartości szczęścia",
     col = "blue", pch = 16)
abline(0, 1, col = "red", lwd = 2)

# Wykres reszt
ggplot(data.frame(resid = residuals(model3), fitted = fitted(model2)),
       aes(x = fitted, y = resid)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red") +
  labs(title = "Reszty vs dopasowane wartości", x = "Dopasowane", y = "Reszty")

model_year <- lm(Life.Ladder ~ Log.GDP.per.capita + Social.support +
                   Healthy.life.expectancy.at.birth + Perceptions.of.corruption +
                    factor(year), data = train)
summary(model_year)

## 
## Call:
## lm(formula = Life.Ladder ~ Log.GDP.per.capita + Social.support + 
##     Healthy.life.expectancy.at.birth + Perceptions.of.corruption + 
##     factor(year), data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.22649 -0.31189  0.02231  0.32743  1.91521 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -0.821926   0.285369  -2.880  0.00408 ** 
## Log.GDP.per.capita                0.226789   0.044129   5.139 3.47e-07 ***
## Social.support                    3.267151   0.238945  13.673  < 2e-16 ***
## Healthy.life.expectancy.at.birth  0.042645   0.006775   6.295 5.09e-10 ***
## Perceptions.of.corruption        -1.604959   0.115291 -13.921  < 2e-16 ***
## factor(year)2014                 -0.021248   0.091286  -0.233  0.81600    
## factor(year)2015                  0.047407   0.091668   0.517  0.60519    
## factor(year)2016                 -0.061083   0.091628  -0.667  0.50519    
## factor(year)2017                  0.029835   0.092085   0.324  0.74602    
## factor(year)2018                  0.106752   0.091941   1.161  0.24595    
## factor(year)2019                  0.096935   0.092418   1.049  0.29456    
## factor(year)2020                  0.089242   0.091499   0.975  0.32969    
## factor(year)2021                  0.055439   0.091926   0.603  0.54663    
## factor(year)2022                  0.011764   0.092790   0.127  0.89914    
## factor(year)2023                  0.108911   0.092071   1.183  0.23720    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5531 on 796 degrees of freedom
## Multiple R-squared:  0.7346, Adjusted R-squared:  0.7299 
## F-statistic: 157.4 on 14 and 796 DF,  p-value: < 2.2e-16

anova(model2, model_year)

## Analysis of Variance Table
## 
## Model 1: Life.Ladder ~ Log.GDP.per.capita + Social.support + Healthy.life.expectancy.at.birth + 
##     Perceptions.of.corruption + Freedom.to.make.life.choices
## Model 2: Life.Ladder ~ Log.GDP.per.capita + Social.support + Healthy.life.expectancy.at.birth + 
##     Perceptions.of.corruption + factor(year)
##   Res.Df    RSS Df Sum of Sq F Pr(>F)
## 1    805 218.43                      
## 2    796 243.50  9   -25.067

AIC(model2, model_year)

##            df      AIC
## model2      7 1251.657
## model_year 16 1357.763

vif(model_year)

##                                      GVIF Df GVIF^(1/(2*Df))
## Log.GDP.per.capita               4.949310  1        2.224705
## Social.support                   2.065306  1        1.437117
## Healthy.life.expectancy.at.birth 3.903559  1        1.975743
## Perceptions.of.corruption        1.200257  1        1.095563
## factor(year)                     1.052572 10        1.002565

pred <- predict(model_year, newdata = test)
mae <- mean(abs(pred - test$`Life.Ladder`))
rmse <- sqrt(mean((pred - test$`Life.Ladder`)^2))
mae

## [1] 0.4325921

rmse

## [1] 0.5405193

error <- test$Life.Ladder - pred
MAE <- mean(abs(error))
MAPE <- mean(abs(error / test$Life.Ladder))*100
RMSE <- sqrt(mean(error^2))
cat("MAE:", MAE, "\nMAPE:", MAPE, "\nRMSE:", RMSE, "\n")

## MAE: 0.4325921 
## MAPE: 8.051326 
## RMSE: 0.5405193

resettest(model_year)

## 
##  RESET test
## 
## data:  model_year
## RESET = 20.643, df1 = 2, df2 = 794, p-value = 1.82e-09

Próba stworzenia modelu uwzględniającego czas. Żaden rok nie jest istotny statystycznie — wszystkie mają p-value znacznie powyżej 0.05. To oznacza, że po uwzględnieniu innych zmiennych (gospodarka, zdrowie itd.), sam rok nie wnosi istotnej zmiany w poziomie szczęścia. Możliwe, że trendy czasowe są już “wychwycone” przez inne zmienne. Anova daje nam jednak informację, że RSS (Residual Sum of Squares) spadło z 545.49 do 532.57, co oznacza, że model 2 lepiej dopasowuje się do danych. Test F sprawdza, czy ta poprawa jest istotna statystycznie. p-value = 0.0004378 → oznacza, że poprawa dopasowania jest istotna na poziomie istotności 0.001. Choć żaden pojedynczy rok w modelu nie był istotny, to rok jako grupa zmiennych kategorycznych wnosi statystycznie istotną informację. Dodanie factor(year) poprawia dopasowanie modelu. Rok wnosi coś istotnego statystycznie, bo nawet niewielka poprawa dopasowania (spadek RSS) przy tak dużej liczbie danych oznacza, że efekt roku jest realny, a nie losowy.

AIC (Akaike Information Criterion) jest również niższy dla model_year. Mimo że model_year ma więcej parametrów (25 vs 7), to i tak wyraźnie lepiej dopasowuje się do danych po uwzględnieniu złożoności. To kolejny dowód (po ANOVA), że dodanie factor(year) ma sens i jest uzasadnione statystycznie.

Mimo że model_year ma lepsze dopasowanie do danych treningowych (niższy AIC, lepszy ANOVA), jego trafność predykcji na danych testowych jest minimalnie gorsza. Różnice są jednak bardzo niewielkie i mogą być przypadkowe (w granicach błędu). Dodanie factor(year) poprawia dopasowanie i ujawnia drobny efekt lat, ale nie poprawia (ani nie pogarsza znacząco) trafności prognozy.

Możemy potraktować tę część projektu jako swego rodzaju eksperyment.

data_2<-dane_finalne
data_2$Region <- countrycode(sourcevar = dane_finalne$Country.name,
                           origin = "country.name",
                           destination = "region")
data_2$RegionCode <- as.numeric(factor(data_2$Region))
levels(factor(data_2$Region))

## [1] "East Asia & Pacific"        "Europe & Central Asia"     
## [3] "Latin America & Caribbean"  "Middle East & North Africa"
## [5] "North America"              "South Asia"                
## [7] "Sub-Saharan Africa"

ggplot(data_2, aes(x = factor(Region), y = Life.Ladder)) +
  geom_boxplot(fill = "skyblue", alpha = 0.7, outlier.color = "red") +
  stat_summary(fun = mean, geom = "point", shape = 20, size = 3, color = "darkblue", 
               position = position_dodge(width = 0.75)) +
  labs(
    title = "Rozkład szczęścia (Life Ladder) w regionach",
    x = "Region",
    y = "Life Ladder"
  ) +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

data_2 %>%
  group_by(year, Region) %>%
  summarise(mean_life = mean(Life.Ladder, na.rm = TRUE)) %>%
  ggplot(aes(x = year, y = mean_life, color = factor(Region))) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  labs(
    title = "Średni poziom szczęścia w regionach na przestrzeni lat",
    x = "Rok",
    y = "Średni Life Ladder",
    color = "Region"
  ) +
  theme_minimal(base_size = 14)

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

anova_region <- aov(Life.Ladder ~ factor(Region), data = data_2)
summary(anova_region)

##                 Df Sum Sq Mean Sq F value Pr(>F)    
## factor(Region)   6  459.8   76.63     122 <2e-16 ***
## Residuals      895  562.0    0.63                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA: p-value < 2e-16 → bardzo silne podstawy do odrzucenia H₀.

To oznacza, że region ma istotny wpływ na poziom szczęścia — przynależność regionalna wyjaśnia dużą część zmienności. Suma kwadratów dla regionu (1205) jest znacznie większa niż dla reszt (1505), co znaczy, że duża część wariancji Life Ladder jest związana z regionem. Regiony bardzo istotnie różnią się między sobą pod względem poziomu szczęścia. To idealne uzasadnienie, żeby uwzględniać region w modelach predykcyjnych lub interpretacyjnych.

9. Podsumowanie (do zmiany)

Model wykazał, że największy wpływ na poziom zadowolenia z życia mają m.in. Log GDP per capita, Social support oraz Freedom to make life choices. Diagnoza reszt wskazuje na pewne odstępstwa od założeń klasycznego modelu liniowego, ale model osiąga stosunkowo niskie błędy (MAE, RMSE), co świadczy o dobrej jakości dopasowania.

W kolejnym kroku możliwa jest rozbudowa analizy o modele panelowe z efektami stałymi lub losowymi, co może poprawić ujęcie efektów krajowych i czasowych.

```{r}}